The onset of tree-like patterns in negative streamers 
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We present the first analytical and numerical studies of the initial stage of the branching process 
based on an interface dynamics streamer model in the fully 3-D case. This model follows from 
fundamental considerations on charge production by impact ionization and balance laws, and leads to 
an equation for the evolution of the interface between ionized and non-ionized regions. We compare 
some experimental patterns with the numerically simulated ones, and give an explicit expression for 
the growth rate of harmonic modes associated with the perturbation of a symmetrically expanding 
discharge. By means of full numerical simulation, the splitting and formation of characteristic 
tree-like patterns of electric discharges is observed and described. 
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It is a well known visible fact that electric discharges 
form tree-like patterns, very much as those in coral reefs 
and snownakes. The study of the branching process 
leading to such pattern is of considerable interest both 
from pure and applied points of view. Many industrial 
techniques, ranging from lasers to chemical processing of 
gases and water purification could be improved provided 
the development of tree-like patterns can be controlled or 
avoided. Although an electric discharge is a very complex 
phenomenon, with radiation and chemistry processes in- 
volved pj[3], the description of its initial stage is simpler. 
A single free electron traveling in a strong, uniform elec- 
tric field ionizes the gaseous molecules around it, gen- 
erating more electrons and starting a chain reaction of 
ionization. The ionized gas creates its own electric field, 
which speeds up the reaction, and a streamer is born. 
The streamers of ionized gas have an inevitable tendency 
to break up at their tips (see figure fl]), followed by the 
creation of the familiar tree-like pattern . 

Early efforts [4-7 were able to identify a minimal 
streamer model with which, after numerical simulations 
under the hypothesis of cylindrical symmetry, an insta- 
bility was observed J5J. Later on, the dispersion relation 
for planar fronts was computed and the existence of an in- 
stability leading to the development of fingers was found 
[9]. Due to the enormous difficulty for performing full 
numerical simulations of the minimal streamer model, 
some simplified descriptions have been attempted in the 
last years (see [10] for a review where various ad hoc 
assumptions are discussed). In any case, the fully 3-D 
case has resisted the attack so far. As an alternative ap- 
proach (the one we follow in this work), the motion and 
propagation of the streamer discharge has recently been 
described by a contour dynamics model first introduced 
in [11] and used to predict with success some experimen- 
tal features of discharges on dielectric surfaces [I2j [13] . 
The contour dynamics model describes the interface sep- 
arating a plasma region from a neutral gas region. For 




FIG. 1. Experimental sequential images of the growth of a 
streamer in a corona discharge at its early stages. This is a 
100 kV point-plane discharge in the vertical direction with 
an air-gap of 30 cm. The contour has been highlighted as a 
visual aid. 



a negative discharge, the separating surface has a net 
charge a and the thickness goes to zero as \[D, with D 
the charge diffusion coefficient. The interface moves with 
a velocity in the normal direction 
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where E^ is the normal component of the electric field 
at the interface when approaching it from outside the 
plasma region, \i e the electron mobility, D e is the elec- 
tron diffusion coefficient, E^ on is a characteristic ioniza- 
tion electric field and k is twice the mean curvature of the 
interface. The parameter Iq is the microscopic ionization 
characteristic length. 

The total negative surface charge density at the inter- 
face changes according to 
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where E u is the electric field at the interface coming from 
inside the plasma, g e is a parameter proportional to the 



resistivity of the electrons in the created plasma and j~ 
is the current contribution at the surface of any source 
inside the plasma. For instance, an insulated wire in- 
side the plasma at Xq, carrying an electric current I(t), 
will create a current density inside the plasma and as 
quasineutrality is fulfilled, we will have for the plasma 
region 



V-j=J(t)(*(x-xo) 



(3) 



and j is obtained solving that equation. Note that at the 
interface there is an electric field discontinuity given by 
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It is convenient to express the model in dimensionless 
units. The physical scales are given by the ionization 
length /q? the characteristic impact ionization field E^, 
and the electron mobility \i e . The velocity scale yields 
U§ = M e Ei, and the time scale To = Iq/Uq. Typical values 
of these quantities for nitrogen at normal conditions are 
l « 2.3 /im, Ei « 200 kV/m, and fi e « 380 cm 2 /Vs. The 
unit for the negative surface density reads do = eoE^/e, 
so for the current density jo = croUo/lo and for the resis- 
tivity £>o = MeW°"o- The diffusion constant unit turns 
out Do = ^o^o- Introducing dimensionless units, the 
model reads 
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(a) Initial configuration 




(b)t = 0.0794 




(c)t = 0.1438 




(d)£ = 0.1769. The encircled region is 
enlarged in figure [3] 
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and £ = D e /Do the dimensionless diffusion coefficient. 
In what follows all the quantities are dimensionless un- 
less otherwise indicated. We have used an adaptive 
boundary element method, developed for general con- 
tour dynamics problems (p~6j [17]) in order to perform 
numerical simulations with equations (pi and (p|. In 
figure [2] we show numerical simulations of the evolu- 
tion of the discharge at four time steps. The plasma is 
assumed ideally conducting, initially charged with inte- 
grated surface charge Q=-25, subject to an external field 
in the vertical direction E = 0.5, and confined inside 
an initially spherical geometry perturbed by r*o (#,(/>) = 
R + (5 (exp(-(cos 2 ((/>) + cos 2 ((9))/c), with c = 0.03 and 
#o = 0.1. We first observe the onset of streamer fingers. 
At time t = 0.17 the streamers develop further instabil- 
ities and split again. Qualitatively the process can be 
described in the following terms: any protuberance that 



FIG. 2. Evolution of the plasma for Q — —25, E 
s — 0.02. Color gradation represents curvature. 



0.5 and 



develops is accompanied by an increase of the charge den- 
sity at its tip. The electrostatic repulsion of charges at 
the tip tends to make the tip expand and the finger grow. 
In opposition to this is the action of the surface tension 
tending to flatten the protuberance and setting up a flux 
of charge from the protuberance out to the sides. How- 
ever, overall, the protuberance becomes amplified. This 
process occurs again and again until a tree-like pattern 
is produced. In figure [3] we depict a detail of this pat- 
tern. Those ideas where anticipated in |5], but due to 
the restriction of 2-D simulations the whole process of 
the branching pattern formation could not be observed. 
In order to be quantitative, we can calculate the growth 
rate of the different modes, both analytically and numer- 
ically. If the initial spherical symmetry is perturbed by a 
small amount, some instabilities will start growing. We 
will study which instability modes are going to prevail 
during the front evolution. 

We consider now a spherically expanding plasma, rep- 
resenting a corona discharge, with Q(t) < so that 




FIG. 3. Detail of the shape of the plasma at t = 0.1769 



E = E+ = 4^)2 < 0. Then, R(t) is given as the 
solution of 

If Q(t) = Q it is easy to check that 
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for the early stages of the discharge and as long as 
R <C ^- . This is in agreement with predictions based 
on continuum streamer models [H],[l5]. If the position 
of the front as well as the charge density are changed by 
a small amount, the perturbed quantities can be param- 
eterized as 
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+ 5V(8,<f>,t), (11) 



where 5 is a small parameter. The angles (9, and <\> are 
the usual spherical coordinates. For convenience we write 
the surface perturbation in terms of spherical harmonics 

as 
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and the surface charge density perturbation as 
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where the coefficients si m (t) and bi m have to be deter- 
mined. Making a standard expansion of the dynamics 
contour model equations (pi and (pi, up to linear terms, 



we get the equations for the particular mode evolution 
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We can get information about the growth of different 
modes by analyzing two special limits. First we study the 
limit of ideal conductivity. It corresponds to g — > 0, and 



hence, from (15), we can conclude that bi m — >• 0. This 



is the case when in the limit of very high conductivity, 
the electric field inside goes to zero (E~ — » 0), as we 
approach to the behavior of a perfect conductor. If we 
consider that Q(t) = Qo is constant or its variation in 
time is small compared with the evolution of the modes 
(which also implies I(t) —> 0), and the same for the radius 
of the front R(t) = ro, we look for a solution si m = 
exp(cjt), cp n = 0, to (14), and get a discrete dispersion 
relation of the form 
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with a maximum at 
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for e <C 1. For a small enough conductivity, g — » oo, 
we find bi m = -E si m (l + 1)/(2Z + 1), and with si m = 
exp(u)i), (fli]) yields 
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for e <C 1. Note that the dispersion relation does not 
depend on m. The finite resistivity cases lay between 
those limits. In fig[4]we have plotted the analytical curves 
given by (16) for different values of e and the results of 



numerical calculations for a perfect conductor. 




and in nature, of the current research in the area of elec- 
tric discharges: bringing the field from a qualitative and 
descriptive era to a quantitative one. 
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FIG. 4. Analytical (in black) and numerical dispersion curves 
for different values of the diffusion coefficient e. The abscissa 
corresponds to the spherical harmonics number, and the or- 
dinate to the growth rate for that mode. 



The dispersion curve allows to predict the expected 
number of branches that will develop. Each branch will 
undergo also a further split and so on propagating to the 
smaller scales. However, it cannot run forever, as there 
is a limitation and the model does not take into account 
the energy radiated, the heat exchange, and the other 
phenomena that will start to play an important role at 
later stages of the discharge. 

The results presented in this work confirm the hypoth- 
esis that at the earlier stages of an electric discharge, 
the main driving forces are diffusion and electrical drift, 
first anticipated in |8]. The expressions obtained for the 



growth rate of the modes, given by (16) and (18) enables 



one to predict the number of forks that one can expect in 
an electric discharge provided the electric field and the 
diffusion constant is known by other means. But, the 
opposite can be worked out: from the numbers of fingers 
observed, we can for example infer the electric field if the 
charge density at the interface and the diffusion constant 
are known. This has been done for the 2-D case [12 . For 
the 3-D case, the density can be obtained from Stark's ef- 
fect measurements, and effective diffusion coefficient may 
be approximately calculated. These results contribute to 
achieve one of the main goals, both in the laboratory 
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